************************************************************************
* Date: October 26, 2017
* David Simon
* Code for matching Schools to Road


/*
1. For each school, calculate shortest distance between each school and  each 
road.    
2. Assign bearing  between road and school
3. merge on Mattis data matched to school with wind direction.  

*/
************************************************************************



set more off
cap log close


	
**************Globals********************
	global home 

	global winddata 
	global roaddata 
	global trafdata 
	global schooldata 
	global output 
	global samples 


log using "$output\allsclaadt25k.log", replace


**********start by Cleaning up the road data**************


* do same thing with US highways
	* the  file make_traffic.do sets up the aadt_latlong.dta file from the
	*FDOT shape files.

use "$trafdata\aadt_latlong.dta", replace


rename (traflat traflong) (latitude longitude)
	gen aadtroad=1
	label var aadtroad "an aadt road"
	
* limit this to roads with an aadt in the top 75 percentile

sum AADT, detail  /*this tells us that the 75 percentile of cars is 18000 */

keep if AADT>=25000

save "$roaddata\aadt25000.dta", replace


***************Create line segments***********************************	
	***Getting prior segment piece's lat/long, if it's from the same roadway	- JAH
foreach var of varlist latitude longitude {
	gen `var'_prior=`var'[_n-1]
	replace `var'_prior=. if `var'==0
	replace `var'_prior=. if `var'_prior==0
	}

	
***getting segments (will be one less than the number of end points (e.g., if you have 3 points, you have 2 line segments) 	
	geodist latitude longitude latitude_prior longitude_prior, generate(segmentdistance) miles
	label var segmentdistance "Road segment distance in miles"
	sum segmentdistance, d
	*hist segmentdistance if segmentdistance>0.1, freq
	
***AND NOW: DIVIDING UP SEGMENTS >0.1: For ordering purposes:
	gen n=_n

	gen segexp=floor(segmentdistance*100)+1  /**If a segment is longer than 0.01 mi,
	round up to the nearest 10th, multiply by 10, and divide it into that number of segments
	Change X in segmentdistance*X to make the minimum distance 1/X miles */
		
	expand segexp
	sort n
	*get a count by segment:
	by n: generate segord=_n

	**Replacing:
	foreach var of varlist latitude longitude {
		gen `var'_original=`var'
		*take the previous latitude, add on the extra fraction of the distance to get to the next value in the data 
		*(or each segment, so multiple by order 1-2-3...)
		gen dist`var' = (`var'-`var'_prior)/segexp if segexp>1
		replace `var'=`var'_prior+segord*(`var'-`var'_prior)/segexp if segexp>1
		*useful for checking: 
		order `var' `var'_original `var'_prior segexp dist`var' segord segmentdistance n
												}	
												
	drop distlatitude distlongitude
	gen id=_n
	
	order latitude longitude segmentdistance aadtroad ROADWAY BEGIN_POST END_POST Shape_Leng id 


	**Chceck 1:  Verify no more long segments:
	***Getting prior segment piece's lat/long, if it's from the same roadway	- JAH
	*foreach var of varlist latitude longitude {
	*	gen `var'_prior_2=`var'[_n-1]
	*	replace `var'_prior_2=. if `var'==0
	*	replace `var'_prior_2=. if `var'_prior==0
	*	}
	***getting segments (will be one less than the number of end points (e.g., if you have 3 points, you have 2 line segments) - JAH	
	*geodist latitude longitude latitude_prior_2 longitude_prior_2, generate(segmentdistance_2) miles
	*	label var segmentdistance "Road segment distance in miles"
	*	sum segmentdistance_2, d
	*	hist segmentdistance_2 if segmentdistance_2>0.1, freq
	***saving this new data as a temporary file:
	
	destring ROADWAY, replace
	
	tempfile aadtrd
	save "`aadtrd'"	

	
********************now create calculate nearest road to school***************
use "$schooldata\NCES_2010.dta", replace

geonear ncessch latcod loncod using "`aadtrd'", neighbors(id latitude longitude) miles //this is where mi_to_nid is made

*merge road data - must set to local drive
rename nid id

merge m:1 id using `aadtrd'
tab _merge


drop if _merge==2
drop _merge



rename latitude road_lat
rename longitude road_lon
rename id road_id



rename ROADWAY roadway
keep road_lat road_lon latcod loncod ncessch mi_to_nid roadway AADT DISTRICT COSITE DESC_FRM DESC_TO AADTFLG





***Make tempfile that will eventually have the list of roads/matches added to it:
clear
tempfile majorrd_matches
gen ncessch=.
save "`majorrd_matches'"

***Identify the number of schools that we'll loop through:
use "$schooldata\NCES_2010.dta", replace
	tempvar all
	gen `all'=_N
	disp `all'


**2. To get the first-fifth nearest road, by school
forv schoollist=1/4298 {  
	forv i=1/5 {
		if `i'==1 {			// First time around, make a temporary file of roads to use for matching; this will be reduced down later
			use "`aadtrd'", clear
			tempfile majorrd_`schoollist'
			save "`majorrd_`schoollist''"
			display "School `schoollist' of 4298"
			}	
		use "$schooldata\NCES_2010.dta", replace // open up NCES list of schools
		
		keep if _n==1 // TO LIMIT TO ONE SCHOOL, EVENTUALLY WE'LL HAVE TO LOOP OVER EACH SCHOOL
		geonear ncessch latcod loncod using "`majorrd_`schoollist''", neighbors(id latitude longitude) miles //this is where mi_to_nid is made

		*merge road data - must set to local drive
		rename nid id

		merge m:1 id using `majorrd'
		tab _merge

		drop if _merge==2
		drop _merge

		rename latitude road_lat`i'
		rename longitude road_lon`i'
		rename id road_id`i'
		rename mi_to_nid mi_to_nid`i'

		rename ROUTE routename`i'
		rename ROADWAY roadway`i'
		keep road_lat`i' road_lon`i' latcod loncod ncessch mi_to_nid`i' routename`i' roadway`i'
		
		local roadtodrop=roadway`i'

		append using "`majorrd_matches'" // append to the list of school-to-road matches
		save "`majorrd_matches'", replace
		
		***NOW, delete the ROADWAY we just matched on and loop back over the next closest one:
		use "`majorrd_`schoollist''", clear
		drop if ROADWAY==`roadtodrop'
		save "`majorrd_`schoollist''", replace
		}
	}

use "`majorrd_matches'", clear

	
**************Now try to calculate initial bearing of school to road  *******************
* check this against this formula:  http://www.movable-type.co.uk/scripts/latlong.html



	g lat1=latcod*_pi/180
	g lat2=road_lat*_pi/180
	g lon1=loncod*_pi/180
	g lon2=road_lon*_pi/180
//=ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1), SIN(lon2-lon1)*COS(lat2)) - from http://www.movable-type.co.uk/scripts/latlong.html
	g theta=atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1)) 
	gen degrees = (180*theta)/_pi
	g bearing=mod((degrees+360), 360)

*reality check: this gives us a bearing between 1 and 360
	sum bearing

	rename bearing angle_degrees



*now I found at least two online calculators to apply the formula to
* http://www.movable-type.co.uk/scripts/latlong.html
* and: https://planetcalc.com/7042/
* and if you want one that does a map with it:  https://www.mobilefish.com/services/distance_calculator/distance_calculator.php

keep if mi_to_nid<=1
	list angle_degrees latcod loncod road_lat road_lon mi_to_nid in 1/10
	
	
	
save "$samples\schooltoaadt25000.dta", replace
